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The quantum magnetic oscillations are studied for planar condensed matter systems with a lin- 
qq \ ear, Dirac-like spectrum of quasiparticle excitations. We derive analytical expressions for magnetic 

I , oscillations (de Haas - van Alphen effect) in the density of states, thermodynamic potential, mag- 

netization, and chemical potential both for zero and finite temperatures, and in the presence of 
i i ' scattering from impurities. We discuss also a possibility of using magnetic oscillations for detection 

, of a gap that may open in the spectrum of quasiparticle excitations due to a nontrivial interaction 

I between them. 
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^3 ! I. INTRODUCTION 

C ; 

Magnetic oscillations (MO) in metals proved to be a powerful tool for investigation of the shape of their Fermi 
i_T surface. The oscillations of magnetization were predicted by Landau in 1930 [2j. The experimental discovery of MO 
came soon: Shubnikov and de Haas (SdH) found oscillations of electrical conductivity in Bi crystals, and later de 
Haas and van Alphen (dHvA) discovered the oscillations of magnetization. Despite more than 70 years old history 
of MO they continue to attract attention of both experimentalists and theorists. The recent experimental studies are 
mostly focused on quasi-two-dimensional (quasi-2D) organic conductors and superconductors (for a review, see e.g. 
£N| '• 0)' The ultimate hope is that better understanding of the 2D and quasi-2D organic systems can also contribute 
OO 1 to the studies of high-temperature superconductors that have a similar layered structure. As suggested recently in 
Ref. Q, dHvA experiment can be used as a probe to detect band- and/or angle-dependent gap amplitudes of type-II 
superconductors. 

It turns out that these experimental advances demand also further development of the MO theory. While MO 
observed in 3D metals are well described by Lifshits-Kosevich (LK) p| (for dHvA effect) and Adams-Holstein @ (for 
SdH effect) theories, there is no commonly accepted and used theory for MO in 2D and quasi-2D materials 0. As 
discussed in Ref. 0, there are two essential features of dHvA effect in 2D case that differ it from 3D case: 
. (1) The sharp saw-tooth like shape of the oscillations seen in the low temperature TCwi and high purity T -C 
regions. Here, lul is the distance between Landau levels and T is the width of Landau's levels due to impurity 
O ■ scattering. 

O ■ (2) Due to the Landau quantization of the 2D kinetic energy in the magnetic field, both the density of electrons, n, 
, and the chemical potential, fi, cannot be fixed simultaneously as the magnetic field is changed. Depending on the 
r physical situation, two extreme limits are usually considered: 

(i) the limit of fixed and field independent chemical potential /i, which can be represented by a grand canonical 
ensemble. This corresponds to the case originally studied by Lifshits and Kosevich 5] in 3D. 

(ii) the limit of fixed density of carriers n, in which a strong dependence of the chemical potential fi on the magnetic 
field that has to be taken into account by considering an equation for /i that leads to the chemical potential oscillations. 
This equation is crucial in 2D case even for the large Fermi energy, tF ^S> wl since despite smallness of the chemical 
potential shift, dfi ~ ljl due to the magnetic field compared to cf, it can change the phase of the magnetization MO 
(in 3D the oscillating part 5/i ~ ujl^l/^f) 1 ^ 2 , so that for ep 3> lol it does not alter the MO with a high accuracy). 
Nevertheless, under certain conditions discussed in Refs. 0Q the oscillations of the chemical potential are small and 
one may still assume that the chemical potential fi coincides with the Fermi energy tp as in the limit (i). 
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On the other hand, in condensed matter systems with a linear, "relativistic" spectrum of quasiparticle excitations, 
the conditions favorable for MO seem to persist to rather small magnetic fields and high temperatures 0. This is 
particularly true for 2D systems with the Dirac quasiparticle spectrum. 

There is a wide variety of planar condensed matter systems that in the low-energy limit have a linear dispersion 
law of the quasiparticle excitations. In particular, the studies of high-temperature cuprate superconductors inspired a 
big interest in the so-called nodal excitations that are the gapless fermion excitations associated with zeros of the gap 
function. Depending on the physical origin of the gap, one can consider rather different physical situations that are in 
general described by different theoretical models. For example, if the gap opens due to the anisotropic electron-hole 
pairing, which is one of the possible states of the electron system with a half-filled band and nested Fermi surface, 
this corresponds to a 2D orbital antiferromagnet (OAF) or staggered flux state [Tol ITU Il2| . 

There is a consensus that the superconductin g st ate in cuprates has a d-wave superconducting energy gap, with 
nodes along the diagonals of the Brillouin zone [LI- A linearization of the Bogolyubov spectrum of quasiparticles 
around four nodes on the Fermi surface leads to another realization of gapless fermion excitations. Although from 
physical point of view the OAF state seems to be more complicated than the d-wave superconducting state, because 
the former is characterized by nonzero local currents violating time-reversal and translational symmetries, the Dirac 
description of OAF state turns out to be simpler since an external electromagnetic field enters into the theory in the 
same way as in QED2+1. 

Carbon-based materials, e.g. pyrolytic graphite and carbon nanotubes present one more example of the planar 
systems with the relativistic dispersion law [Til ITU Il6| . From an experimental point of view these materials are 
probably the most promising because the quality of available samples already allows one to observe SdH effect 0, 
and to study quantum effects that are due to a high magnetic field 0] (see Ref. [2(j for a review). 

The relativistic theories in an external magnetic field have been the subject of research in quantum field theory for 
many years (for reviews see Ref. The extraction of MO from t he g eneral expressions presents a rather subtle 

and interesting problem that was investigated in Refs. |22l 123. l24l l25l l26| . However these studies of the MO in QED 
are mostly concentrated on the field theoretical aspects. 

The purpose of the present paper is to make a systematic study of the MO in QED2+1 devoting special attention 
to the link with planar condensed matter systems and to quantities that are particularly important for condensed 
matter theory, e.g. density of states (DOS), thermodynamic potential and magnetization. To make the theory more 
realistic we also include into the model the effect of scattering from impurities, by considering Landau levels with field 
and temperature independent width. 

We begin by presenting in Sec. |n] three 2D models mentioned above that can be studied using the same effective 
relativistic Lagrangian written in Sec. Ill Dl in the presence of an external magnetic field. The Green's function 
necessary for subsequent calculations is introduced in Sec. IIIII and the DOS oscillations are studied in Sec. IIVI both 
with and without scattering from impurities. The general representation for the thermodynamic potential in terms of 
the DOS and its link to the corresponding nonrelativistic potential are considered in Sec.[V] In Sec. I VII we derive the 
expression for the thermodynamic potential at T = in the absence of impurities with explicitly extracted MO (the 
calculational details and an alternative representation for the thermodynamic potential in terms of the generalized 
zeta function are given in Appendices 1X1 and |B"|) . In Sec. I VI Al and at the end of Sec lVIliTl we discuss a possibility of 
detecting a gap that may open in the spectrum of one of the systems discussed above. The analytical expression for 
the magnetization at T = in the absence of impurities written in terms of Bernoulli polynomials is given in Sec. lVIII 
In Sec. I Villi we generalize the results obtained in the previous sections to the presence of impurities (Sees. I Vlll Al 
and IVIIIBl with the calculational details presented in Appendix and in Sec. IVIII CI for nonzero temperature. In 
Conclusions, Sec. IIXI we give a concise summary of the obtained results. 

II. MODEL 

There are many planar condensed matter models that in the low energy sector can be reduced to QED2+1 form 
and here we briefly describe some of them underlining the main assumptions leading to the Dirac-like form of the 
effective Hamiltonians. 
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A. Hamiltonian of d-density wave state 



Due to d x 2_ y 2 momentum dependence of the gap in the OAF states it is also called rf-densi ty w ave (DDW) state 
[27j • The mean-field phenomenological Hamiltonian describing this state can be written as 0, [2^ 

ff DDW =E / 7^£(k)k% 3 -/ 1 i-#)4,(k), (2.1) 
Jrbz (2tt)^ L J 



<T=U 



where the spinors 



Xa(k) = ( c Jk ( +Q) )' X f CT (k) = (4(k) 4(k+Q)), (2.2) 



are composed from creation and annihilation operators cj.(k), c CT (k) for momentum k and spin cr, the single particle 
energy is e(k) = —2t(cosk x a + cos k y a) with t being the hopping parameter, /1 is the chemical potential, D(k) = 
-^■(cosfc^a — coskyd) is the c£-density wave gap and Q = (n/a,Tr/a) is the wave vector at which the density-wave 
ordering takes place. The integral is over the reduced Brillouin zone RBZ and Ci are Pauli matrices. Throughout the 
paper h = c = fc^ = f units are chosen, unless stated explicitly otherwise. 

This Hamiltonian describes the excitations with the spectrum EQs) = —fx ± \J e 2, (k) + D 2 (k) . Linearizing the 
spectrum about the four nodes Ni = (±7r/2a, ±7r/2a) with i = 1, . ..,4 at half-filling (fi = 0) one obtains E(k) = 

— ^ ± \Jv F k x + Vpky, where the Fermi velocity calculated for half-filling vf — |<9e(k)/<9k|k=N| = 2\/2ta, the DDW 

gap velocity V£> = |9Z3(k)/9k|i c= N| = -^D a, and the momenta k x and k y are given in the local nodal coordinate 

system (see Fig. 2 of Ref. [H). 

Using this linearized spectrum and inserting the vector potential in a way that preserves the charge conservation 
for the original Hamiltonian (|2.1() , one can arrive 29] at the following Dirac Lagrangian describing the quasiparticles 
in the vicinity of one of the nodes, e.g. N^ = i = (7r/2a, 7r/2a) 

£ = xUxWD v xUx), x=(t,r), (2.3) 

where Xai x ) = X l Ji x ) <J i is the Dirac conjugated spinor, i labels the node N^, the covariant derivatives D v are 

fhdt — ieAo(x), v = 0, 

v F (hd x -ilA x {x)), u=l, (2.4) 
va {Kdy - i- c Mix)) , v = 2, 

and the 7 matrices are 

j v = (a 1 ,-ia 2 ,ia 3 ), {f, 7"} = ^hg^, gT = diag(l, -1, -1), fi,v = 0,1,2. (2.5) 

The chemical potential fi can be introduced in the Lagrangian i|2.3f) by choosing Aq = /i/e and the sum over the spin 
components in Eq. I|2.3|l can be regarded as an additional flavor index. 

In 2+1 dimensions, there are two inequivalent representations of Dirac algebra (see e.g. 

7° = cr 3 , j 1 =i<Ji 1 l 2 = i<J2, (2.6a) 
7 = -cr 3j 7 1 = -icri, 7 2 = -icr 2 , (2.6b) 

which correspond to right and left handed coordinate systems. As one can check, our 7-matrices (|2.5ll correspond to 
the representation l|2.6b|l if one makes a unitary transformation x 1 = Ux 1 , 7 M = U'j^U^ 1 with U = (l/2)(/2 — io\ — 
i(T2 + 103) Since the physical properties of the system depend only on the algebra that these matrices obey, one can 
directly work with the more commonly used representation l|2.6|l instead of (|2.5|) . 

Dealing with l|2.3|l one should not forget that all physical quantities involve the summation over the 4 nodes 
present in the original Hamiltonian l|2.1|l . In fact, only two neighboring nodes, e.g. N^ = i = (ir/2a, ir/2a) and 
Nj = 2 = (n/2a, — n/2a) are nonequivalcnt in RBZ and the corresponding local nodal coordinate systems are related 
to each other by parity transformation (k x ,k v ) — > (k x ,—k y ) with the simultaneous interchange vf <-> vd- Thus the 
sum over nonequivalent nodes can be taken into account by doubling spinors 

T. = (|) (2-7) 
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and using the reducible representation of 7-matrices, 7^ = (03, ia%, ia 2 ) ® 03: 



with the following Lagrangian density 



£ = T (T (x) l7 l/ ^T ff (x). (2.9) 



B. Hamiltonian of d-wave superconducting state 

In contrast to the hypothetical DDW state, d-wave superconductivity (<iSC) is observed in cuprate superconductors 
[l3| . It can be described by the Bogolyubov - de Gennes Hamiltonian 



J dW(i,r)r 3 [e (-iV - r 3 ^A(r)) - /,] *(i,r) 

/" d 2 n J d 2 r 2 [A t (t,r 1 ;r 2 )* t (t,r 1 )r_*(t,r 2 ) + * t (t,r 1 )T + *(t,r 2 )A(t,r 1 ;r 2 ) 



(2.10) 



where 'J and ll^ = (V>j,V ; l) are the standard Nambu spinors, e(p) is the single particle energy, t± = [t\ ± ir 2 )/2 
and A(t, ri; r 2 ) is the bilocal gap operator the Hermitian conjugate of which includes the transpose in the functional 
sense, i.e. At(i,n;r 2 ) = [A(i, r 2 ; ri)]*. The nontrivial dependence of A(i, ri;r 2 ) on the relative coordinate ri — r 2 
is supposed to describe d-wave pairing state. Further simplification of the pairing term of Eq. I|2.10l) is possible if 
one assumes that the amplitude of the pairing term is a constant with respect to t and the center of mass coordinate 
R = (ri +r 2 )/2. For cuprates this assumption may well be justified even above the critical temperature, T c . Below T c 
in the vortex state this assumption corresponds to neglecting vortex core contributions and considering pure "phase 
vortices". Since we deal with the amplitude of the bilocal complex field A(t, 1*1, r 2 ), it is natural to consider also its 
phase 9(t, R). There are, however, some subtleties in writing the continuum version of the pairing term of 12.1l)fl in 
an explicitly gauge i nvariant form that can be solved by choosing an appropriate form of the differential operator for 
A (see Refs. |3lt l32j). Since we look for the low-energy quasiparticle excitations near four gap nodes, we write the 
linearized Hamiltonian for one of these nodes as (see e.g. |28l l31j ') 

H dS c = ^Hx) [iv F T 3 (d x - i-T 3 A x ^ + iv A T + e ie{x ^ 2 d y e l0(x ^ 2 + iv^e-^^dye^ 9 ^' 2 ^ V(x), (2.11) 

where vf is the Fermi velocity, «a = |<9A(k)/<9k|ic=N| is the gap velocity defined by the slope of the superconducting 
gap A(k). Both velocities are calculated in the nodal points on the true Fermi surface of the system, but not at 
half-filling as for the DDW case. The linearization puts restrictions on the domain of validity of the Hamiltonian 
(|2.11ll which, however, remains rather wide |33| . 

Finally, introducing the Dirac conjugated spinor ff = 4 ,l t T2 we ca n obtain from the Hamiltonian (|2.11|) a Lagrangian 
similar to Eq. (|2.3|) . but with the vector potential that couples only with the vf term. Again making, if necessary, 
a unitary transformation, one can combine spinors rf originating from two opposite nodes into one four-component 
spinor. 

Thus in the absence of the field the low lying quasiparticle excitations in the dSC are described by the relativistic 
dispersion law 



E(k) = ± y Jv F k 2 +vikl (2.12) 

where the momenta k x and ky are given in the local nodal coordinate system. Experimental values of vf and for 
cuprates can be found in Ref. |34| . for example, the value of vf ~ 2.5 x 10 5 m/s and the anisotropy of the Dirac cone 
vf/v\ varies between 10 and 20 depending on doping and compound. 

There are two important differences between DDW and <iSC Hamiltonians. Making a gauge transformation that 
removes the phase from the last term of Eq. I|2.11|) , we observe that quasiparticles in superconductor do not couple 
simply to the vector potential A corresponding to the magnetic field B = Vx A, but to the supercurrent V9—(2e/hc)A, 
where 9 is the phase of the superconducting order parameter. This difference between DDW and dSG cases leads to a 
conclusion that in an external magnetic field the nodal quasiparticles of DDW state form Landau levels |3f| , while in 
dSC state Landau levels are strongly mixed j3|J. Nevertheless we presented dSC Hamiltonian ]2.11|l here due to its 
practical importance and the fact that the simplified physical picture based on the Landau levels |37l l38|| may become 
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relevant for higher energies |39j and can be regarded as the first approximation to a more complicated case of the 
vortex state in d-wave superconductors. 

The second difference comes from the fact that for DDW state chemical potential /i is explicitly present in the 
Lagrangian (|2.3|) . while for dSG state it was absorbed in the definition of the Fermi velocity vf on the Fermi surface. 
The origin of this difference can be traced back to the different structures of the Hamiltonians 1|2.1[) and (|2.1(J|1 . One 
can also say that in dSC state the chemical potential of nodal quasiparticles is zero, so that when the applied field is 
changed the corresponding Landau levels (even if they were formed) cannot cross the chemical potential and produce 
dHv A effect H3- 



C. Layered graphite 

The semimetallic e nerg y band structure of a single graphene sheet has the conduction and valence 7r bands with 
the energy dispersion [4(| 




E(k x ,k y ) = ± £ yi + 4cos^^cos^+4cos 2 ^, (2.13) 

where t « 3eV , a = \fiacc — 2.46A is the lattice constant of two dimensional graphite, acc is the distance between 
two nearest carbon atoms. These two bands touch each other and cross the Fermi level in six K points located at 
the corners of the hexagonal 2D Brillouin zone, but only two of them, for example, K = (27r/a)(l/-\/3, 1/3) and 
K' = (27r/a)(0, 2/3) are inequivalent due to the periodicity of the Brillouin zone. The low-energy excitations can be 
studied by taking the continuum limit a — > at any two independent K points labelled as j — 1,2. They have a linear 
dispersion Ek — ±i>f& with vp — (V3/2)ta ~ 9.7 x 10 5 m/s. These excitations can be formally described by a pair of 
two-component (Weyl) spinors ipj a , which are composed of the Bloch states residing on the two different sublattices 
of the bi-particle hexagonal lattice of the graphene sheet. The Hamiltonian describing these excitations, for example, 
in the point K coincides with the free Dirac one 0, 0] 

H = E /(^^(^l^^+fW^^k): ( 2 - 14 ) 

where the momentum k = (k x ,k y ) is already given in a local coordinate system associated with a chosen K point, 
"01(7 = *Pi a J and j ' 1 ' 2 = (cr 3 , i<j2, — ivi)- This representation of 7-matrices can be mapped in the representation 
(|2.6a[l by using a unitary transformation with U — (l/\/2)(I + iuz)- Since the local coordinate system for the 
point K' is related to the system associated with the point K by a parity transformation, these two spinors can be 
again combined in one four-component Dirac spinor ^ a = {ip\ a , ip2a)- The number of spin components N has to be 
regarded as an adjustable parameter and N — 2 corresponds to the physical case. Finally, the Lagrangian density of 
noninteracting quasiparticles reads 

£0 = V v F ^(t, r) f *r°(ft + iM) _ n ig x _ h 2g \ <^ (t r)j (2.15) 

where ~ ^P^ and 4x4 7-matrices are either ((73, 103, —iai) <g) 03 |4l| or can be taken from the unitary equivalent 
representation Q2.8[l [bH |42|. Since the terms with d x ^ y in Eq. (|2.15|) originate from the usual kinetic term of the 
tight-binding Hamiltonian, vector potential A can be inserted in the Lagrangian (|2.15|l using a minimal coupling 
prescription. Finally we note that three-dimensional version of the Dirac Hamiltonian (|2.14() was used in Ref. |9j to 
describe an unusual magnetoresistance in the two doped chalcogenides Ag2+aSe and Ag2+,5Te. 



D. Model relativistic Lagrangian 

As we have seen, many models of planar condensed matter systems result in the Dirac type form of the effective 
low energy theory. Thus as a starting point of the present paper we choose the following Lagrangian 



£ = Qiii'fDp - A)^, M = 0,l,2, 



(2.16) 
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= is the Dirac conjugated spinor and the vector potential 
for the external magnetic field B perpendicular to the plane is taken in the symmetric gauge 



where = <9 M — ieA^ is the covariant derivative, ip. 



x I B B 



(2.17) 



The nonzero chemical potential /x will also be taken into account by choosing Aq = (i/e, in the energy- momentum space 
this corresponds to a shifting u> — > u> + fx. We assume that the fermions carry an additional flavor index i = 1, . . . , N 
which can be used to calculate the sum over equivalent nodes in the case of dSC case and to sum over the spin 
components in the cases of DDW and graphite. In (|2.16|) we have already set the velocities c = vp = vd = v& = 1, 
so that we consider the "isotropic Dirac cone" . When it is needed they can be restored according to the prescription 
discussed in Refs. |2flll4^. The Dirac 7-matrices are taken in the reducible four-component representation and they 
obey the Clifford (Dirac) algebra, 7^7" + 7^7^ = 21^ iV . The explicit form of the 7-matrices is given in Eq. (|2.8(1 and 
the symmetry properties of the Lagrangian with reducible four-dimensional representation are discussed, for example, 
in [3(J, 0i ■ We have also included a mass (gap) term A in the Lagrangian l|2.16|l . The physical origin of this gap 
depends on the underlying system we consider. For example, it is well known that an external magnetic field is a 
strong catalyst in generating such a gap for Dirac fermions (the phenomenon of magnetic catalysis) |44| . Usually the 
opening of the gap marks an important transition which occurs in the system. In particular, in the case of pyrolytic 
graphite a poor screening of the Coulomb interaction may lead to excitonic instability resulting in the opening of the 
gap in the electronic spectrum and manifesting itself through the onset of an insulating charge density wave (see e.g. 
|4lL l42j). There are, however, many obstacles in experimental detection of such a gap, so that in the context of the 
MO studies we consider in Sees. IVI Al and IVIII CI a possibility for its detection using dHvA and SdH effects. 

A full theory of MO should also include the Zeeman interaction term which leads to a spin factor in the LK formula. 
Here, however, we neglect this term motivated by the fact that for the relativistic spectrum the MO may become 
observable for the relatively low fields when the spin splitting is still small. If necessary, the Zeeman term can be 
added explicitly both to the original Hamiltonian (|2.1() and the Lagrangian (|2.1(i|) . 



III. SPECTRAL FUNCTION OF DIRAC QUASIPARTICLES IN AN EXTERNAL FIELD 

The Green's function of Dirac fermions described by the Lagrangian H2.16|l in an external field given by the vector 
potential H2.17|l reads 



S(x — y) = exp I ie J A\dz ] S(x — y), 
v 



(3.1) 



where S is the translation invariant part S of the Green's function. Its derivation using the Schwinger proper-time 
method and decomposition over Landau level poles has been discussed in many papers (see e.g. Refs. |29|, |44j, |45|L 
so that here we begin with the spectral function associated with the translationary invariant part S of the Green's 
function 



A{w,p) 



1 

27ri 



S A (u-iQ,p) -S R (u + iO,p) 



(3.2) 



where the retarded, S R , and advanced, S , Greens' functions are written in the energy-momentum representation. 
This spectral function decomposed over Landau levels reads 01 



/ 2 \ ^ 

A(«, P )=e«p(-iU5>ir 



( 7 °M n + A)/ 1 (p) + / 2 (p) 



2M n 



5(lu - M n ) + 



( 7 °M n - A)A(p)-/ 2 (p) 
2M n 



5(lo + M n 



where M n = VA 2 + 2eBn and 



(3.3) 



A(P) 



eB 



P+L n -! 2 



eB 



/ 2 (p) = 4(p 1 7 1 +p 2 7 2 )LLi (2 



(3.4) 



with P± = (1 ± sgn(e_B)i 7 1 7 2 )/2 being projectors and L n , L\ Laguerre's polynomials (L 1 _ 1 = 0). In what follows, for 
convenience, we take eB > 0. The chemical potential //, as mentioned above, has to be taken into account via the 
shift lo — * lo + /I. 
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In contrast to the nonrelativistic Landau levels with the energies E n — ^^-(n + 5) (here m is the effective mass 
of the carriers), in the relativistic problem with zero gap A the energies are E n = ^j ehv ^ v ^ B2n, n — 0,1,..., 
where we restored all parameters (V2 = vp,vo or va depending on the model we consider) to show explicitly the 
differences. For the nonrelativistic problem the distance between Landau levels coincides with the cyclotron frequency, 
lu c = ehB/(mc), so that w c [K] ~ 1.35i3[Tesla]r7i e /m, where m e is the electron mass, while for the relativistic problem 
the corresponding energy scale, characterizing the distance between Landau levels, is 



uo L = J hVFV22eB [K] = 4.206 x 10-4 S^m/sy^Tesla], (3.5) 
V c y vf 

where vp is given in m/s. For example, choosing m e /m = 1 we estimate that u> c ~ IK • _B[Tesla], but this value can 
be increased by using metals such as Bi with a large ratio m e /m. Making the estimate for the relativistic case with 
a rather small vp — 2 x 10 5 m/s and vp/v2 = 20, which are typical values of the parameters for the DDW model, we 
obtain u>l ~ 18A" • -J B [Tcsla] which shows that in the systems with the linear dispersion law the quantum condition 
favorable for MO persists to rather high temperatures and small fields 0- Moreover, repeating the estimate for the 
Fermi velocity in graphite, vp = 9.7 x 10 5 m/s, we obtain even larger ujl ~ 4001^ ■ \J B [Tesla] . We note that since the 
Zeeman term has the same magnitude as u> c , it is indeed small comparing to u>l and can be safely neglected. 

In order to consider the MO for a more realistic case, one should introduce the effect of quasiparticle scattering 
that results in a Dingle factor in the expression for the amplitude of MO. In general, this can be done by considering 
dressed fermion propagators that include the self-energy S(w) due to the scattering from impurities. Up to now the 
problem of scattering from impurities in the presence of a magnetic field does not have yet a satisfactory solution. 
Therefore, here we choose the case of constant width T = T(u) = 0) = —ImE R ((J = 0) = l/(2r), r being a mean free 
time of quasiparticles, so that the 5-like quasiparticle peaks corresponding to the Landau levels in Eq. I|3.2J) acquire 
a Lorentzian shape: 

Such a broadening of Landau levels with a constant T was found to be a rather good approximation valid in not 
very strong magnetic fields 0, 0] . Definitely, the treatment of disorder in the presence of the magnetic field in such 
a simplified manner should be considered as only the first step until further progress in this problem is achieved 
(in connection with this, see, Ref. |47j). The technical advantage of the approximation T = const is that one can 
firstly consider dHvA oscillations with (5-like Landau levels at T = and only afterwards introduce the effect of level 
broadening due to the finite temperature and quasiparticle scattering convoluting the final results with the appropriate 
distribution functions (see, Refs. 0, EH an d Sec. 0). However, this simplification is only valid if all Landau levels 
have the same width. 



IV. OSCILLATIONS OF DENSITY OF STATES 



We begin with the calculation of the quasiparticle density of states (DOS) Do(e) in the absence of scattering (r = 0), 
which is expressed in terms of the spectral function (|3.3|l as 



Do(e) 



d 2 p 
(2tt) 2 



tr [y A(e,p) 



(4.1) 



where the domain of integration B is chosen to preserve the volume of the original Brillouin zone. We also assume 
that the summation over spin states that are identical in the absence of Zeeman term is included in the definition of 
trace in D(e). This, as mentioned above, can be done by taking an appropriate value of N. 

Evaluating the trace and expanding the limits of integration over momenta to 00, we get the DOS as the sum of 
delta functions of Landau's level energies: 



D (e) = 



NeB 
2tt 



S(e -A) + 5(e + A) + 2j2 W e - M„) + *(e + M n )) 



n=l 



(4.2) 



The fact that the original Brillouin zone has a finite volume can be taken into account later using the finite limits of 
integration if necessary and/or by implying that the DOS D(e) is multiplied by the factor #(A — |e|), where A is the 
bandwidth. 



8 



The broadening of Landau levels due to the scattering is taken into account according to the prescription given by 
Eq. (|3.6|l . It is easy to see that this prescription corresponds to the convolution of the density of states Do(e) with 
the Lorentz distribution P r (uj) = T/[tt(lli 2 + T 2 )] 



oo oo 

D(e)= J duPc(u)-e)D (u), J du>Pr(u) = l. 



(4.3) 



In fact, instead of the Lorentz distribution any other normalized probability distribution can be used [l| in order to 
introduce phcnomenologically the energy level width. In what follows we shall use the Lorentz distribution which 
allows analytical calculations. Physically, the Lorentz distribution corresponds to an extremely strong disorder since 
all moments of this distribution diverge. The quasiparticle DOS at the Fermi surface can be obtained by evaluating 
Eq. i|4.3[l for e = [i. Eq. (|4.2|) can be also rewritten as 



A>(e) 



iVeB|e| 
NeB 



2n Sgn(£) de 



S(e 2 -A 2 ) + 2J2 <S(e 2 - A 2 - 2eBn) 
d 



(e 2 - A 2 ) + 2 - A 2 - 2eBn) 



n=l 



(4.4) 



where is the step function. Note that NeB /{2k) is the density of the Landau levels. Using the Poisson summation 
formula 



oo oo 
1 oo „ oc „ 

-F(0) + Y J F(n)= / F(x)dx + 2ReY, / F(x)e 2nikx dx, 

71=1 { k=l{ 



(4.5) 



we find the sum over the Landau levels 



n=l 



2eBn) = 6{e z 



1 e 2 



2eB 



E 

fc=i 



sin ( irk 



eB 



(4.6) 



Substituting now the last expression in Ea. (|4.4l) we obtain the final expression for the DOS with zero broadening that 
can be written in three equivalent forms: 



a*) = ^«6»w| yd 2 - a 2 ) 



2 A 2 o ^ 1 • /7rfce 2 -A 2 

e - A + 2eB > — sin — i '- 

^ nk \ eB 



k=l 



2eB 



■ tan cot 



2 _ a - 
~eB 
n(e 2 -A 2 ] 



2eB 



e 2 - A 2 - 2eBB! ' ' 



2eB 



The summation in Eq. (|4.7bl) was done using the formula 

OO . / N 

Esin (Trnx) , 
— = tan" 1 

n 



sin(7ra;) 
1 — cos(7rx) 



(4.7a) 

(4.7b) 
(4.7c) 

(4.8) 



and Eq. (|4.7c|) is written using the Bernoulli polynomials B n (x) periodically continued beyond the interval [0, 1]: 

B n (x) = £ J^r cos (2^ - ^) , n>2, 0<x<l; n = 1, < x < 1. (4.9) 



fc=i 



2 / 



In fact, all the Bernoulli polynomials we dealt with in this paper depend on the fractional part mod[x] of their 
argument x, i.e., B n (mod[x\) (here mod[x] is a shorthand notation for x modulo 1, i.e., mod[x] = x — [x] where [x] is 
the largest integer satisfying [x] < x). In what follows, for brevity we omit the sign mod[] and write simply B n (x). 
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A. DOS in the presence of scattering 



Now we turn to calculating DOS in presence of the impurity scattering rate T using convolution (|4.3|l . Since -Do(e) 
is an even function of e, the expression for D(e) can be written as follows 



N 

D(e) = — Im { (e + iT) 



A /eB 



dx 



dx 



X + Z IT J (x + z) 



tan 



( cot y) 



eB 



A 2 - (e + iT) 2 
eB 



(4.10a) 



(4.10b) 



where we used Eq. (I4.7b(l for -Do(e), represented Pr as Pr(e±uj) — — Im(e±ui + iT)^ 1 and integrated by parts to obtain 
the second term in square brackets in Eq. I|4.10a|) . We also put an explicit cutoff A associated with the bandwidth in 
the first integral describing non-oscillating part of DOS. 
It is easy to calculate that 



Im 



A /eB 



(e + iT) 



dx 

x + z 



= Tin- 



A 2 



i/(e 2 - A 2 - T 2 ) 2 + 4e 2 r 2 



+ e tan 1 



2eT 



A 2 



(4.11) 



The second integral in Eq. (|4.10|) denoted as / can be calculated if we take into account that the function 
tan -1 (cot(7ra/2)) is periodic with the period x = 2 so that we can write 



1=1 



dx 



2fc+2 



7T J (x + z) 2 



■E 



tan ^cot T j=-^ 



fc=0 



dx 



tan 1 ( cot — 



2 k 



{x + zf 



TlX 



k=0 



dx 



(x + 2n + z)'- 



■ tan 



i 

1 (cot y) = \J dx{1 " 2x)C(2 ' x + z/2) ' 



(4.12) 



where £(s, x) is the generalized zeta-function and we used also that tan 1 (cot(7rx/2)) = (7r/2)(l — x) for < x < 2. 
The last integral is evaluated using the properties of ^-function 



J dx((2,x + z) = ~, J 



dx xC(2, x + z) = *(1 + z) - In z, 



(4.13) 



so that we finally get 



/=-,i-)+ln|-i. 



Thus, we express the DOS in terms of the digamma function ip: 

'A 2 - (e + ^^) 2 



D(e) 



or, in equivalent form 



N 

72 



A 2 

r In Im 

2eB 



(e + iT) U 



N A 2 
De = ilrin 



2eB 7T 2 de 



-Im 



lnT 



2 eJ B 



( e + l r) 2 



2eB 



eB 



-In 

2 



A 2 - (e + iT) 2 
A 2 - (e + T 1 ) 2 



2eB 



(4.14) 



(4.15) 



(4.16) 



It is evident that the DOS oscillations are contained in the V'-function when the real part of its argument becomes 
negative. They can be extracted in explicit form using the relationship 



1p(—z) = 1p( z ) H 1" 7TCOt(7Tz). 

Z 



(4.17) 
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Hence we get 



A 2 - (e + iTf 



Re?/> 



2eB 
c 



A 2 - {e + iT)' 



2eB 



(e 2 -A 2 -r 2 ) + 0(A 2 +r 2 -e 2 )) 



2 -A 2 -T 2 |-2ier\ 



2eB 



i sgn(e 2 - A 2 



r 2 )Im7/> ( ' 



2eB 



(e 2 - A 2 - T 2 ) 



2eB 



e 2 - A 2 - T 2 + 2ieT 



IT COt 7T- 



2eS 



so that Eq. I|4.15[) is rewritten as 



In 



A 2 



2eS 

+ esgn(e 2 — A 2 — T 2 



ReV' 



ft 



A 2 -r 2 |-2zer 



eSle 2 - A 2 - T 2 



Imi/j 



2 eJ B 



A 2 -r 2 | -2ier 



2eS 



(e 2 - A 2 -T 2 ) 2 + 4e 2 r 2 
2eBeT 

+ ( e 2 - A 2 - r 2 ) 2 + 4e 2 r 2 



Tr9(e 2 



r 



2 ^esinh(27rer/eS) - rsin[7r(e 2 - A 2 - T 2 )/eB] 
cosh(2ireT/eB) - cos[7r(e 2 - A 2 - T 2 )/eB] 



(4.18) 



(4.19) 



For small T the main contribution comes from the last term in Ea. (|4.19|l . thus DOS in the presence of T can be 
represented in the form similar to Eq. (|4.7all : 



D(e) = f sgn(e)^ {6(e 2 - A 2 - T 2 ) [ £ 2 - A 2 - T 2 
2tt 



d 
'd~c 



2eBj2 



k=l 



1 7rfc(e 2 - A 2 

— — sin — 

7tk eB 



■ exp 



27rfc|e|r 



(4.20) 



Using now the formula 



E—e kt sin fca; = tan 
h 



sin .t 



t > 0, 



k=l 



cos a; 



(4.21) 



one can check that the last term in Eq. I|4.19(l is recovered. It is evident that oscillating part of DOS is contained in 
the sum over k in Eq. (|4.20|) . 



B. DOS in limiting cases 

Since the final expressions for DOS, Ea. I|4. 15|) and especially Ea. (|4.19|l are rather lengthy, it is useful to consider a 
few simple limiting cases and compare them with the known results. First of all, it is easy to obtain from Eq. (|4.15(l 
that in the limit of zero field B = the DOS becomes 



m = S 



rin- 



A 



v / r 2 + (e-A) 2 



Tin- 



A 



x /r 2 + ( e + A)^ 



7T 

— + tan 



2er 



(4.22) 



and for A = 0, after restoring the prefactor 1/(vfva) with the velocities vf and va, it reduces to the DOS derived 
in Ref. HJ. 

As follows from Eq. (|4.19|l . the DOS at zero energy, but in the finite field is given by 



AT 

D(0) = — 

IT" 



In- 



A 2 



r 2 



- * 



2eB 



In 



A 2 



eB 



2eB 



r 2 



(4.23) 



The first term of Eq. (|4.23(l is nothing else but the zero energy DOS (|4.22(1 in the absence of the magnetic field. The 
behavior of the DOS (|4.23|l can be now studied in various asymptotical regimes. For example, for A = and r — > 
we find 



D(0) = — 



eB A 2 
- +rln 2^ 



(4.24) 
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i.e., DOS is enhanced in the presence of the magnetic field if T 2 -C eB. In the opposite limit, eB <C T 2 we have 



D(0) 



NT 



hi 



A 2 



r 2 T 3r 4 



On the other hand, for A ^ and r^Owe obtain 

D(0) 



NT 
T 2 " 



' A 2 ( A 2 
In — — - * 



2eB 



2eB 



eB_ 
A 2 



(4.25) 



(4.26) 



so that nonzero gap regularizes T — > divergence, which is present in Eq. Q4.24[) . 



V. REPRESENTATION FOR THERMODYNAMIC POTENTIAL 



A. General expressions 



All thermodynamic quantities such as the magnetization M and the number of electrons N can be found from the 
grand thermodynamic potential £1, which in the relativistic case (see Refs. |25j,|44J) can be written using the DOS as 
follows 



oo 

tl(T,n) = -T J deD(e)\n ^2 cosh 



2T 



(5.1) 



with the DOS D(e) given by Eq. (4.3). We assume everywhere that the volume (area) of the system is unity, so 
that Eq. I|5.1|l corresponds to the thermodynamic potential per unit volume. It is convenient to separate the vacuum 
contribution at T = 0, fi = in the thermodynamic potential. For that we write 



n(T,/t) = -T / deD(e) In ( + ) (0(e) + 0(-e)) 



deD(e)(e - ^)sgn(e) — T / deD(e) 



In 1 + e~ 



In 1 



(5.2) 



or, using the evenness of the function D(e): 

oo 

0(T,/x) = - J deeD(e) -T j d< Di< ) 
o 



In 1 



9(e) + In ( 



1 + e~ 



(5.3) 



The first (divergent) term in the last expression is the vacuum energy while the second one (convergent) is due to 
contributions of real quasiparticle excitations. 
At zero temperature, T = 0, we thus have 



fi(0,A*) = n(0,A* = 0) + J tfeD(e)(e-||i|)=fio(0) + f2o(A0. 
o 

It is easy to see that the density of states is related to the thermodynamic potential at zero temperature 

d 2 Q{T = 0,/i) 



dy? 



(5.4) 



(5.5) 



As was already mentioned, the constant Landau level width approximation appears to be very useful simplification, 
viz. instead of calculating the potential (|5.1(l with f ^ 0, one can start from the potential fir, 



n T (fi) = -T I deD (e)ln 2 cosh 



2T 



(5.6) 
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where Do(e) is the DOS with zero level broadening (T = 0). Then the potential fl at finite T is obtained by convoluting 
Q,T{fJ-) with the distribution function Pr(u>). The statement that the level broadening is equivalent to considering the 
distribution of the chemical potentials was proven in Ref. |48j and we adopt here its derivation for the relativistic 
thermodynamic potential Ea. (|5.1|) . Indeed, if the levels are broadened, then the density of states is obtained as a 
convolution of Do(e) with the probability distribution Pr(e) of energies e: 



oo oo 

Q(T,n) = -T J dwdeP T (uJ - e) D (w)ln ^2cosh^^=-T J dude Pr(ui)D (e) In ^2 cosh 



6 — LJ — /l 

2T 



(5.7) 

oo oo V" / 

= —T J dujdeP r (uj- n)D (e)\n ( 2cosh -^r^ = ^ (iw^rfw-^rH, 

— oo — oo 

where the potential ^It(^) is given by Ea. (|5.6(l . If several damping effects occur together, the corresponding convo- 
lutions have to be carried out successively (the order is not essential). 

The effect of finite temperature can also be included in this scheme by choosing the distribution function Pt(z), 
which describes the temperature line broadening, equal to the negative derivative of the Fermi function, i.e. 

Pt(z) = = Ij—. (5.8) 

W dz 4Tcosh 2 ^ V ; 

We now show that the thermodynamic potential Eq. (|5.6[) can be obtained as the convolution of fto(n) with the 
distribution function Pt given by Eq. (|5.8|l 

oo 

O r (/i)= J dwP T (u- ju)fio(w). (5-9) 

— OO 

Since the vacuum contribution does not change under averaging over thermal and T distributions, it is enough to 
prove this fact for the finite part ^It(^)- Performing integration by parts we obtain: 

oo oo e 

Q T (H)=-T J deD{e) In (l + e 1 ^) 0(e) + In (l + e^) 0(-e) =- J def F (e) J dxsgn(x)D(x) 7 (5.10) 

— OO — OO 

where /f(c) is the relativistic thermal distribution function 

/F(e) = ^L + JfciL, (5.U, 

1 + e t 1 + e t 

which is also used in semiconductor physics. Integrating by parts for the second time, we obtain 

oo 

OrM - / dnF %~^ W(e)de, (5.12) 



where 



€ X 

W(e) = J dxsgn(x) J dy sgn(y)D(y). (5.13) 
o o 

Since 

-^ = ^TT--*W. ^0, (5.14) 
oe 4T cosh ^ 

we find that 

W{n) = -n (M), (5-15) 



13 



hence 

oo 

Or(M)= / de(- dnFi g~^ ^) ()„(,). 

— OO 

The common effect of level broadening due to both temperature and damping effects can be written as successive 
convolutions 

oo 

fl(T,fi)= I dLo'dujP T {J - fi)P T {uj-uj')n {oj). (5.17) 



Hence for calculating thermodynamic quantities we need to know only the thermodynamic potential Qo(u>) at zero 
temperature and zero width of levels that is directly expressed via the DOS -Do(e), Eq. Ij4.7|l . Thus the knowledge of 
zero T (and zero T if the Landau levels have the same width) DOS is completely sufficient to write down the finite T 
thermodynamic potential. 

B. Link with nonrelativistic thermodynamic potential and equation for chemical potential 

We start from the expression for the nonrelativistic thermodynamic potential (see, e.g. Ref. Q) 

oo 

r* N R,(T» = -T J deV{e)ln(l + e li ^y (5.18) 

— oo 

expressed in terms of the DOS T>(e) for original nonrelativistic tight-binding Hamiltonians considered in Sec.|n] Due 
to the complexity of the tight-binding spectrum, this DOS cannot be found explicitly for B ^ and the only property 
of the DOS we need is that it is an even function of e. 

The derivative of the thermodynamic potential (|5.18() with respect to the chemical potential /i, 



dn NR (T,y) 



deV(e)n F (e- fi) = n, (5.19) 



determines the density of carriers n in nonrelativistic many-body theory as a function of T, B and /i. On the other 
hand one can consider Eq. I|5.19[) as an equation for /i as a function of T, B and n which is typical for studies in a 
canonical ensemble. In Eq. (|5.19(l np(uj) = l/[exp(cj/T) + 1] is the usual Fermi function. At T = 0,/i = Eq. (|5.19(l 
gives the density of particles for a half-filled band. 



DC 

J deD{e)9(-e), (5.20) 



so that the deviation of the particle density from hq is due to a finite temperature and nonzero chemical potential. 

It is convenient to redefine the thermodynamic potential in order to have an explicit proportionality of (i to the 
density of free carriers. Thus we introduce 

oo 

Q'(T,n) =O NR (T,M)+/m = - J deV(e) [r\n(l + e^ - ^(-e)] . (5.21) 

— OO 

Inserting 1 = 0(e) + 9(— e) before the logarithm, the last expression is rewritten in the form 

oo oo 

n'(T,fi) = - J deeV(e)-T J deV(e) In (l + e 1 ^) 9(e) + In (l + e^) 0(-e) , (5.22) 

-oo 
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where the first term gives the energy due to the half-filled zone. As is seen, fi'(T, p.) acquires the form of a relativistic 
thermodynamic potential (compare with Eq. I|5.3[) ). Its derivative with respect to p 



^'f^ = J Me - /*)0(e) - [1 - n F (e - p)}0(-e) 

— CO 

1 f°° 6 

deV(e) tanh ■ 



(5.23) 



2 ./_ „ v ' 2T 



determines the charge density of carriers or carrier imbalance, p (p = n — no = n + — rt_ , where n + and n_ are the 
densities of "electrons" and "holes, respectively) as a function of T and B at fixed /x. Note that now p — at p = 
and there is a symmetry with respect to the transformation p — > —p and p — ► —p. In the present paper, however, we 
will mostly use the canonical ensemble interpretation of Eq. (|5.23|) . viz. as the equation for p as a function of T and 
B at fixed p. 

VI. ZERO TEMPERATURE THERMODYNAMIC POTENTIAL 

The zero temperature thermodynamic potential is given by Eq. I|5.4|l , and for calculation of its vacuum contribution 
^o(O) we refer to Appendix [X] The thermodynamic potential with vacuum energy f2o(0) subtracted is given by the 
second term of Eq. (|5.4|l 

M 

n (M) = / deD (e)(e-\n\), (6.1) 
o 

where the DOS for T = case is given by Eq. 1)4. 7fl . In what follows we choose for definiteness p > 0. Calculating 
Clo(p) (see Appendix IB1 for the detail) we obtain that it can be represented as a sum of regular and oscillating terms, 

CIq(p) = VL reg (p) + tt osc (p), (6.2) 

with Q reg (/x) expressed in terms of generalized zeta function 



Vl reg (p) = - — 0{p - A) 



N 

and 



\p{p 2 -<Stf)-AeB-{2eBf/^(-\> 1 



A 2 



(6.3) 



N(eBf/ 2 ^ 1 

^ osc (At) = 2^ ^ ~ A ) (7rfc) 3 / 2 [^^^ cos(Trfcw) + J 2 (nkv) sin(7rfcxx>)] , (6.4) 

where the monotonic functions J1.2 are defined by Eq. I|B5(1 . Writing the last expression for Q osc (p) we introduced 
new variables 

a 2 - A 2 a 2 

«;= — , "=-5. (6.5) 

eis eB 



For small fields, e£> << p 2 , we can use the asymptotic expansion (|B7Jl for Ji, J2 and represent r^ osc in terms of 
periodically continued Bernoulli polynomials defined in Eq. I|4.9|l : 

In particular, keeping the first few terms in the expansion we have 



N(eB 



,2 



n osc (p) ~ — i-0(p - A) 

Z7T/1 



„ /"w\ en (w 
B > ( 2 ) + 3^ 3 (2 



eB /w\ (eB) 2 /w 
' i?4 I — 



(6.7) 



3a* 2 a V 2 7 ' 4/i 4 * V 2 . 

where the explicit expressions for the Bernoulli polynomials B\ - B4 are 

B 1 (x)=x-l/2, B 2 {x) =x 2 -x + l/6, B 3 (x) = x* - 3x 2 )2 + x/2, B A (x) = x 4 - 2x 3 + x 2 - 1/30. (6.8) 

The terms oscillating with a magnetic field (Eqs. I|6.6|l . lE7Jl) represent small corrections to the nonoscillating part of 
the thermodynamic potential fEq. H6.3fl ) at eB <C p 2 ■ Nevertheless, since Q osc contains fast oscillating functions of the 
variable (p 2 — A 2 )/2eB, after differentiating it over the magnetic field, the magnetization, M = —d£l/dB, acquires a 
large oscillating part (see Sec. IVIII and Pigs. ^ HI)- 
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A. The period of oscillations and its dependence on A for fixed p 

As it follows from the expressions i|6.4|) , i|6.6[l , the frequency of oscillations of the thermodynamic potential is equal 
to w/2 = [p? — A 2 )/2eB. Therefore their period SB satisfies the relationship 

where lul is defined in Eq. (|3.5[) . Using the estimate of ojl for the DDW model made below Eq. (|3.5[) . we obtain that 
for A = 0, p = lOOA' and B = lTesla, the period of oscillation 5B = 0.032Tesla, which is clearly within the detectable 
range of dHvA oscillation experiments and close to the estimate SB w 0.05Tcsla made in Ref. 35] . 

Repeating this estimate with B = lTesla for graphite with p = 200A" and A = 0, we get SB = 4Tesla which is 
very close to the observed value for period of the SdH oscillations in this material ■ However, so large value of SB 
clearly violates the condition SB <C B implied to obtain Eq. I|6.9[) . and for the case of graphite the period of MO is 
given by a more accurate expression 



1 1 



eB; eB 



p 2 — A 2 A ' 



1 (6.10) 



where Bi and are the values of the magnetic field in two adjacent minima. In particular, for graphite again using 
for an estimate p = 200K and A = we obtain \l/Bi — l/Bi + i \ = STeskiT 1 . Writing the last equality in Eq. HO.lOfl . 
we introduced the area of an extremal cross section of the Fermi surface, A = 7r(/i 2 — A 2 ) as it is usually done. Thus, 
contrary to conventional systems with quadratic dispersion law and the period of MO ~ 1/p, the period of MO in 
the systems with linear dispersion (for A = 0) is ~ 1/p 2 as seen in Fig. ^ 

Making above estimates we have assumed everywhere that A = 0. However, as it is clear from Eqs. 1)6. 9f) and 
(|6.10(l . the opening of the gap A increases the period of MO. In particular, for graphite even a rather small value 
of A, e.g. A < 10~ 2 ^, would produce a sizeable ~ 0.04Tesla change in the period of MO. The effect of the gap 
opening on the oscillations of the magnetization is shown in Fig. [21 This figure was obtained using the values of the 
parameters typical for graphite (see the caption). As one can see, the gap opening produces an observable change in 
the period of MO. Thus this method can be useful to test the realization of magnetic catalysis phenomenon |44j in 
graphite (see Ref. H2)i when the external magnetic field can induce the opening of an insulating gap A(B) in the 
relativistic spectrum of quasiparticle excitations. The crucial condition for the observation of the SB increase caused 
by the gap opening is that the chemical potential p itself does not change 50] as the gap A opens. Formally this 
condition corresponds to considering fixed p in the grand canonical ensemble. However, if the experimental setup 
corresponds to the fixed carrier imbalance p, one can see from Eqs. (18.31) and l|8.4l) below that the entire difference 
p 2 — A 2 — 2ir\p\/N adjusts to the strength of the applied field, so that the period of MO remains unchanged, making 
the gap detection impossible. Nevertheless, as we discuss at the end of Sec. IVIII CI there is still a possibility of the 
gap detection if simultaneously with a frequency the amplitude of dHvA oscillations is measured. 

As pointed out in Ref. (35|, the estimate of SB is not sufficient for the detectability of dHvA effect. In addition we 
should calculate the magnetization and estimate the magnitude of its oscillations. 

VII. MAGNETIZATION 

The magnetization M in the direction perpendicular to the plane should be calculated in the canonical ensemble, 
thus we need the Helmholtz free energy which is related to the thermodynamic potential as 

F(p, B) = ilQiip, B),B) + M (p, B)p. (7.1) 

Then 

dF(p,B) dn(p,B) } 

M p (p,B) = ^— = ——\^ {pB) = M fl (p(p,B),B), (7.2) 

where is the magnetization obtained in the grand canonical ensemble at constant /i: 

M M B ) = Qg— | M =const> ( 7 - 3 ) 

and we used that dfl(p, B)/dp = —p. Formally, the magnetization calculated in the canonical ensemble has the same 
form as the magnetization calculated in the grand canonical ensemble, but we have to take into account oscillations 
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of p in the system with fixed p. In the present paper we restrict ourselves in what follows by studying MO in the 
limit lo c <C e_F, so that for fixed p, as shown in Sec. IVIII Al one can neglect the oscillations of p(B). 

The thermodynamic potential consists of the sum of three terms, the vacuum energy IJA4JI . regular (|6.3|) and 
oscillating (|6.4|l parts. We begin with the calculation of the part of the magnetization which is due to the vacuum 
energy Ea. HA4() . 

As one could see from Eq. I|5.22|) . this energy corresponds to the energy of the half-filled zone. Thus this case is 
also interesting from physical point of view [IJ. Considering firstly the limit A = 0, we obtain from Eq. I|A5(1 that 
the magnetization 

M(A = )M = 0) = _ 3JV C(3/2) e 3/ 2 ^ (? 4) 

4V27T 2 



(note the diamagnetic character of the vacuum contribution in magnetization) . The ~ VB dependence of the magne- 
tization at fx = implies that the susceptibility x = OM/dB oc _B -1 / 2 diverges at zero field raising a question about 
the stability of a homogeneous state described by effective Lagrangian l|2.16[l . As discussed in Refs. |l2l l35| in the 
case of the OAF (DDW) this instability disappears if a coupling between layers is included, or if finite temperature 
and/or chemical potential is included. The instability is also removed in the presence of the nonzero gap A. Indeed, 
using the asymptotics for zeta function 

(M + o) ^__L + _J_ (7 . 5) 



for a^> 1, we obtain from Eq. (|A4|) that 

Np p B 

M(A, M = 0) = -^^-, eS«A 2 . (7.6) 

Thus the opening of the gap removes the singularity at B = in the susceptibility. 

Now we turn to the case p > A when in addition to the vacuum contribution two more terms appear in the 
thermodynamic potential, the regular Q reg (p) and oscillating fl osc (p) parts, that at T — are given by Eq. (|6.3|) and 
Eq. i|6.4|) . As one can easily see, for p > A field dependent terms in Eq. (|6.3|l coincide with the corresponding terms 
in the vacuum energy (|A4|) up to a sign, so that their total contribution to the magnetization is zero. Therefore, only 
the oscillating term of the thermodynamic potential given by Eq. (|6.4|l contributes to the magnetization. 

In the small field limit, eB <C p? differentiating the potential l|6.6|l we obtain 

Nep0(p-A) ~ r(n+l/2) r /«a ,ws un (2eB\ n+1 

M) ~ 2^ („ + !). F+ 2 { 2 J Bn+1 I 2 ) 2 \ {-jT ) ( 7J ) 

where we used the relation dB n (z) / dz — nB n ~\{z). Ea. l|7.7JI is the generalization to the relativistic spectrum of the 
formula for the magnetization obtained by Shoenberg [l| . For example, keeping first few terms in the expansion 17. 7|) 
we have 



M(A = 0,/i) 



Nep 
2tt 



/v\ 3 eB /u\ le 2 B 2 /v\ „ f e 3 B 
* -2^(2) -2^ 3 (2) 



(7. 



where the explicit expressions for the Bernoulli polynomials are given in Eq. (|6.8|l . We recall that the Bernoulli 
polynomials used here depend on the fractional part of their argument so that magnetization is oscillating function 
with the period SB given by Eq. 1)6.9(1 and weakly varying amplitude. The dependence of the magnetization on the 
inverse field B^ 1 is presented in Figs. and [21 computed on the basis of Eq. I|7.7|l . The dependence M(B^ r ) has a 
clear saw-tooth like shape showing that for the considered range of the parameters, it is indeed convenient to represent 
the oscillations in terms of the Bernoulli polynomials. 



VIII. INFLUENCE OF IMPURITIES AND TEMPERATURE ON MAGNETIC OSCILLATIONS 

A. Thermodynamic potential Cl(p) and equation for p in the presence of scattering 

To obtain the equation for the chemical potential it is convenient to use the DOS expressed as in Eq. H4.16|l . which 
allows one to write down the thermodynamic potential l|6.1|l , including the effect of impurity scattering, as 



o (m) 



^rinii 

2tt 2 



NeB 



2eB 



de Im 



lnT 



(e + iTf 



2eB 



1 / A 2 -(e + zr) ^ 

2 n l 2eB 



(8.1) 
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Correspondingly, Eq. I|5.23[l for the carrier imbalance at T = and finite T takes the form 



dQ {p) Nn_. A 2 JVeB T 
= — — 1 In — lm 

dp 7T 2 2eB it 2 



A 2 -(p + iT) 2 \ 1 fA 2 -(p + iT) 



lnT — - + -In 



2e5 / 2 V 2eB 



•2) 



For small r -C p, A, one can derive from Eq. 18.211 a more simple expression 



, 2 - a 2 I 2eB tan-i sin[2^ 2 -A 2 )/(2e£)] , 2 

7j- exp[2^r/(ei?)]-cos[27r( At 2 -A 2 )/(2 e i?)]' ^ >A ' ^ 

where we have introduced the Fermi energy ep, counted from the edge of the gap A, by means of the relation 
(N/2n)(e F — A 2 ) = p. Note that Eq. (|8.3() can also be derived directly from the relationship p = JT deDo(e) and 
Eq. I4.2()j) using the sum I4.21|l . 

In general, as mentioned in the Introduction, both fixed p and fixed p cases are possible. For example, fixing p 
in Eq. (|8.2() one can study the oscillations of the carrier imbalance p as a function of the field B. Eq. (I8.3fl can be 
considered as the relativistic analog of the corresponding relation between p and ep in nonrelativistic case JjJ. One can 
deduce formal "nonrelativistic" limit from Eq. (|8.3[) if one defines the nonrelativistic Fermi energy ep(nr) = np/NA 
and corresponding "nonrelativistic" chemical potential p — p nr + A and consider p nr << A. This nonrelativistic limit 
is, in fact, irrelevant for MO considered in the present paper, because "relativistic" theory studied here emerges as an 
effective theory for nonrelativistic models described in Sec. [HI so that it is meaningless to consider its nonrelativistic 
limit. 

On the other hand, for the magnetization the case of a fixed carrier imbalance p is more natural, and therefore 
oscillations of the chemical potential have to be taken into account by solving Eq. (|8.3|l for p(B). Generally, this 
nonlinear equation cannot be solved analytically without further approximations, but for ep 3> eB the second term in 
Eq. (|8.3|) represents a small correction to the constant value of the chemical potential so that we obtain by iterations 
an approximate solution 

„2 = e 2 _ 2eB , sm[27r(e| - A 2 )/(2ei?)] 

7r exp[2^ F r/(eB)]-cos[2 7 r(e 2 ,- A 2 )/(2eS)]" [ 1 

It is plotted in Fig.|21for the model parameters typical for graphite. As expected, these oscillations have saw-tooth like 
shape and the sharpness of the teeth depends on the value of T. One can also anticipate that the relative amplitude 
of the oscillations decreases as the Fermi energy ep increases and the field B decreases. This is indeed seen in Fig. [3] 
Yet in 2D the chemical potential oscillations are rather strong even for large ep, due to the particular form of Eq. 1)8. 3|) 
that turns out to be different from its 3D counterpart. 

In spite of the smallness of the second term in Eq. (|8.4I) , it is important to take into account the difference of p and 
ep in the arguments of oscillating functions (in the magnetization, for example) in 2D, because it alters the phase 
of oscillations. The difference can be neglected if the additional condition 2-KpT > eB is valid (see discussion about 
the importance of the chemical potential oscillations in Ref. @). This problem, however, is beyond the scope of the 
present work. 



B. Thermodynamic potential Q oac (p) and magnetization M osc in the presence of scattering 

We return now to the oscillating part of the thermodynamic potential l|6.1(l including the effect of impurity scattering. 
The oscillating part can be written (see Appendix lOl for details) as 

^ osc (^ r) = [C ' 9(p 2 -A 2 - T 2 ) M 3/2 [ J i i^ky, nk-i) cos(irkw) + J 2 (irkv, nkj) sin(Trfcw)] , (8.5) 



H (7^) 3/2 



where the functions of two arguments Ji(p,r) and J2(p, f) that generalize the functions Ji and J 2 defined in Eq. (|B5p . 
arc defined in Eq. (|C4|) and the variables w, u and 7 are 

2_ A 2_ r 2 A 2 +r 2 r 2 

w= , u= — , 7 = ~5- ( 8 - 6 

eB eB eB 

Note that the variable w (and also u) that was defined in Sec. I VII bv Eq. (|6.5() is redefined in the presence of scattering 
in accordance with Eq. 1)8. As before the frequency of MO is given by to/2, so that it would diminish in the presence 
of impurities if we considered fixed p case. This decrease of the MO frequency for fixed p is specific for the present 
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relativistic model, while in the nonrelativistic case the MO frequency remains intact even in the presence of impurities. 
The origin of the difference between relativistic and nonrelativistic cases can be traced back to the difference in zero 
field DOS for these cases. 

For r, \feB~ « \i keeping the leading term in asymptotic expansions I|C8J| . (|C9|) (this is the first term in Eq. @Sj 
only), we obtain from Eq. I|8.5|l (omitting theta function) 

N(eB) 2 ^ cos(7r/cw) / (jT 



It is easy to check that for T = Eq. (|8.7|l reduces to the expression that follows directly from Eq. I|tj.4[l if one 
substitutes the first term of the expansion l|B7(l , J\ (p) i=s 1/ v /p. This comparison with Eq. (|6.4J) allows to identify in 
Eq. I|8.7|l the Dingle factor 

Rn(k, At) = cxp (-Zirk^J , (8.8) 

which describes reduction of the amplitude of k-th harmonic due to the electron scattering. As is seen, the finite 
amount of impurities leads to smoothing of the oscillations. The factor l|8.8|l is a generalization of the Dingle harmonic 
damping 

R D R {k) oc exp (-2nk^-J , (8.9) 

for the nonrelativistic case, usually written in terms of the nonrelativistic cyclotron frequency uj c . Comparing Eqs. i|8.8|) 
and 1)8. 9|) . one can notice two facts. First one was already discussed in Sec. lIIII viz. in the relativistic case Rd can be 
larger than it would be for the nonrelativistic case for the same field for /i < mvpV2 = uj\/(2uj c ), so that the condition 
for the observation of MO in the relativistic case is less strict. In particular, for graphite this condition remains valid 
for fj, < m e v F 6 x 10 4 K. Secondly, further increase of /i makes MO unobservable, contrary to the nonrelativistic 
case, when Rd does not depends on the chemical potential //. 

Calculating the magnetization (|7.3H and keeping only dominant terms, we obtain from Eq. (|8.7|) 

A^ 2 -A 2 -r 2 )^sin(7rfc W ) _ 2nkM r 



2nu ^—^ nk 

k=1 (8.10) 
N(fi 2 - A 2 - T 2 ) _! sin[27r( M 2 - A 2 -r 2 )/(2e^)] 

~~ 2^ m exp[27TA(r/(eS)] - cos[2tt(a< 2 - A 2 - T 2 )/(2eB)\ ' 

where in the second equality we used Eq. I|4.21|l . We will show now that the oscillating part of the magnetization 
M osc is directly related to the oscillating part of the chemical potential. 

Comparing the result Ea. (|8.1U|l for M osc with the solution for chemical potential l|8.4|) . one can notice that the 
oscillating part of the magnetization is directly proportional to the oscillating part of /j, 2 : 

M„„ = ^A^VW = (8-11) 

(this is obvious if one neglects small terms T 2 in the trigonometric functions in Eq. I|8.10|l 'l. This proportionality 
can be in fact seen when Fig. [3] (plotted for finite T) is compared with Figs. ^ and [21 (plotted for T = 0). Note that 



unlike nonrelativistic case |jj the magnetization is proportional to n 2 sc , not /i osc - If the ratios e 2 F /eB and 2iry,T/eB 
are small, the oscillating part of the chemical potential under trigonometric functions in Eq. 1)8. 10|) leads to additional 
frequencies in the Fourier spectrum of the magnetization. In its turn this changes the harmonics amplitudes compared 
with the LK like formula Eq. I|8.10|) . In the non-relativistic case this was demonstrated in Ref. 0, • 

C. Calculation of Q, osc (n) for T/0 

Now we consider the effect of finite temperature on the oscillations of the thermodynamic potential that can be 
calculated from the zero temperature potential (|8.5() using the convolution property Eq. (|5.9() . For that it is convenient 
to rewrite the expression (|8.5|l in the form 



fe=l v ' 



Vi{t + 1) 



.12) 
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where we used that the functions Ji(p, r) and J2(p, r) can be represented as — Im and Re parts of the same function 



y/nJi(p,r) = -Im 



dte - pt -r/t 



Vi(t + i) 
u 

and rotated the integration contour to the imaginary axis. Hence 



\pn J 2 (p, r) = Re 



dt e-vi-rlt 
Vt{t + i) ' 



(8.13) 



N(eB) 



3/2 



27T 3 / 2 



E 



i 



J M0 3/: 



-Im 



dt e i7rfc T/* 
\/*(* + 1) 



dc 



4Tcosh 2 



(8.14) 



x 9(e 2 



T 2 ) exp 



itnke 2 
eB 



ink- 



- A 2 - r 2 



At low temperatures, T — ► 0, after making a shift e — > e + /i and changing the variable of integration e — > 2Te, the 
integral over e gives (neglecting T 2 term in the exponent) 



dc 



cosh 2 e 



e -w(*+!) 4T ^ = 2 



cosh 2 e 



/47>(t + l)7rfc 



4ir 2 kTfi(t + 1) 

2ir 2 kT f*(t+l) 
eB 



eB sinh 



(8.15) 



where we used the formula (3.982.1) from Ref. [52j. We obtain 
fi osc (T,/x,r) : 

V" 

OO 



N{eBfl*T» » 2 a 1 



x Im 



e 4 e 



-inkw 



sinh a* 2 * WD 



(8.16) 



Nonoscillating factor in the integrand has a maximum at t = 0, thus we can take sinh evaluated at t — while the 
remaining integral over t is evaluated by means of Eqs. £U), G3- We finally get (omitting again theta function) 



n osc (T,n,T) =NeBTj2 



cos(nkw) 



1 



Trfc S i n h 27r2fc7 > 

k=l bUm eB 



exp 



JV(eB) 2 v cos(fb) 



.17) 



where we introduced the temperature amplitude factor 

2Tr 2 kTn/(eB) 



R T (k,n) 



sinh^ 



(8.18) 



which is a "relativistic" equivalent of the temperature reduction factor of the famous LK formula. Clearly, since 
Rr{k,lj) — > 1 for T — > 0, Eq. I|8.17|l reduces to Eq. (|8.7|) and for finite T the thermal broadening causes a reduction 
of the MO amplitude with respect to T = case. This becomes more evident if we consider the limit B — > 0, so that 
the function 1/ sinh[27r 2 fcT/i/(ei?)] ~ exp[— 2-K 2 kT [ij{eS)\. Comparison of this exponent with Eq. (|8.8|l allows one to 
introduce the Dingle temperature 



T D = -, 



(8.19) 



so that a reduction of amplitude due to quasiparticles scattering from impurities can be interpreted as leading to an 
effective rise of a temperature from true temperature T to T + Td, i.e., as if the system could not be cooled below 
a minimal temperature Tjj. Note that while the Dingle factors l|8.8|l and (|8.9|l are different, the value of Tjj itself is 
the same for relativistic and nonrelativistic cases. Comparing Eq. I|8.17|l with the corresponding expression (15) for 
flosc in 2D from Ref. (this expression contains also a spin factor R s which is not considered in our work), one can 
notice that relativistic and nonrelativistic cases differ only by the factor (~l) fe that appears due to the presence of 
the zero energy, Eq — uj c /2 for the nonrelativistic Landau levels. 
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It is seen from Eq. (|8.17|) that the conditions for observation of dHvA oscillations can be formulated as an inequality 
for the strength of magnetic field: 

M ~ A > eB > max{27r 2 l>, 2n 2 T D p}. (8.20) 

The first inequality is necessary in order to have at least one oscillation. It implies physically that the field must not 
be so high that all fermions are in the lowest Landau level, i.e. the filling factor has to be bigger than one. The second 
inequality guarantees that the amplitude of oscillations is not exponentially suppressed, and for that the magnetic 
field must be strong enough to make spacing between Landau levels bigger or of the order of the thickness of the 
thermal layer times the Fermi energy. 

The dependence of the magnetization (|7.3(l on B^ 1 is obtained by differentiation of Eq. I|8.17|) and the results of 
numerical computations are shown in Figs. 01 and [S] They illustrate the above mentioned facts that the MO are 
smoothed as the temperature and/or the chemical potential rise. Obviously the increase of T gives the same result. 
Moreover, one can clearly see the damping of dHvA amplitude as the field B decreases. It is known that the damping 
of dHvA oscillation amplitude as function of several parameters is commonly used to extract system parameters such 
as Dingle temperature and the effective electron mass. However, these methods of analysis depend on the applicability 
of the LK theory. Assuming that the relativistic generalization of the LK theory considered here is valid at least for 
lol "C £f, we suggest that, using Eqs. (|8.10|) and (|8.17|l . one can extract both the value of p from the damping and 
the difference, p 2 — A 2 from the frequency of dHvA oscillations. The knowledge of these two parameters allows to 
obtain the value of the gap A even for the fixed carrier imbalance p. 



IX. CONCLUSIONS 



In this paper we have studied magnetic oscillations in thermodynamic quantities such as DOS, thermodynamic 
potential and magnetization in planar systems with the relativistic Dirac-like spectra for quasiparticle excitations. 
The attention was mainly paid to the regime, where the calculations in canonical (fixed p) and grand canonical (fixed 
p) ensembles give equivalent results. Our main results can be summarized as follows. 

(1) We have obtained analytical expressions that describe MO in the DOS given by Eqs. (|4.7() . (|4.19() for zero width 
r = (no impurities) and Eq. I|4.20|l for finite width T for Landau's levels. 

(2) For r = 0, we have expressed the MO in thermodynamic potential, Eq. (|6.6|l and magnetization, Eq. I|7.7|l as 
a series in the periodically continued Bernoulli polynomials. This representation turns out to be useful for 2D case 
when the oscillations have saw-tooth like shape. In the limit u>l <C p it is sufficient to consider only first few terms 
of the series. 

(3) For finite impurity scattering rate T and temperature we have obtained the thermodynamic potential (|8.17l) and 
represented it in terms of the Dingle factor Rp, Eq. I|8.8I) . and the temperature amplitude factor Rt, Eq. H8.18fl . as 
it is usually done in Lifshits-Kosevich theory. The spin factor, R s , can be calculated similarly if the Zeeman term is 
included in the model. 

(4) For finite T, we have derived also the equation (|8.2|) for the chemical potential. Its solution for fixed carrier 
imbalance exhibits oscillations of p as a function of the magnetic field. It is shown that the oscillating part of the 
magnetization is in fact directly proportional to the oscillating part of the chemical potential (more precisely, to 

(m 2 W)- 

(5) On the basis of obtained formulas, we have discussed the possibility to detect a gap that may open in the spectrum 
of the Dirac-like quasiparticle excitations due to a nontrivial interaction between them. 

One of the possible and necessary extensions of our results would be to study magnetic oscillations for the case when 
the consideration within the canonical ensemble is crucial, i.e. in the regime where the oscillations of the chemical 
potential described by Eq. I|5.23|l at fixed p (cf) cannot be neglected. The necessity of such an extension emerges 
from the fact that the oscillations of p become important for 2TrpT < eB which is exactly the condition (|8.20() for 
having large amplitude of the MO. 

Another important extension would be to consider MO analytically in the transport quantities such as electrical 
and thermal conductivities. These oscillations were in fact seen in the numerical results presented in Ref. |29j, but 
analytical results can be useful for comparison with the experiments 0, 0] . While for dHvA effect the condition 
p = const is more natural, it is plausible that SdH effect can be measured under condition p = const. This, as we 
have discussed in the paper, can be used for an experimental observation of the gap A, especially if the gap depends 
on the applied field, A = A(B). 

All above-mentioned problems can be treated analytically due to the fact that the broadening of Landau levels 
has a Lorentzian shape and the impurity scattering rate T (the width of this distribution) is assumed to be field and 
temperature independent. In fact, both these assumptions can be questioned, in particular for 2D systems with a 
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linear dispersion. For example, the validity of the first assumption for 2D systems is now discussed in the literature 
(see Refs. in Q). While the second assumption may well be valid in the low field regime, it is necessary to investigate 
its domain of validity and probably to consider also the dependence T(B) to access the high-field regime. 
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APPENDIX A: CALCULATION OF VACUUM ENERGY Q (fi = 0,B) 

As mentioned above, the vacuum contribution does not change under averaging over thermal and T distributions, 
so it is sufficient to calculate it at T = and in the absence of impurities. We calculate the vacuum term using the 
density of states in the form of Eq. I|4.2|l 



DC 

n (n = 0,B) = - J deeD {e) 



NeB 
2tt 



A + 2 ^2 M n 



n=l 



(Al) 



The sum over the Landau levels is divergent and to calculate it we use the representation 
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to write 
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-leBt 



n r dt 
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eBt coth(eBt), 



(A3) 



where we introduced the ultraviolet cutoff A as the bandwidth. The integral can be evaluated through the generalized 
zeta function 



n (p = 0,B) 



N 
2^ 



AA 2 , N , /a / 1 A 2 



O 



(A4) 



where the term ~ A 3 which does not depend on A and B is omitted. The last expression coincides with the 
vacuum energy computed in the second paper in Ref. 0]. Setting in Eq. (|A4|) A = and using the identity 
C(-l/2, 1) = C(-l/2) = -(1/4tt)C(3/2) we obtain 



n (^0, Bl A = 0) = ^^(eB)^. 



(A5) 



To compare Eq. I|A5(1 with Eq. (15) of Ref. 01 an d Eq. (10) of Ref. [3j|, one should restore the prefactor 1/(vfvd) 
with the velocities vf, vd and express them via the parameters t and Do as done below Eq. H2.1|l . Note that the 
ground state energy increases with applied magnetic field, indicating diamagnetism. 

APPENDIX B: CALCULATION OF n (/x) 

Substituting Eq. (|4.7a|l in Eq. I|6.1(l and integrating by parts we obtain 
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where 



N 



and 



fio 0*) = -g^O" - - A) 2 ^ + 2A), 



^ > M = - -—6{n A) ^ — / de sin ■ 
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N{eBf' 2 y, i r dxsm{Tikx) 
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Writing the last expression for r4 2 V) we introduced new variables w defined in Eq. H6.5(l and u — A 2 /(eE). 
To extract explicitly the oscillations in the variable w we rewrite the integral in Eq. <|B3|) as 
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where v is defined in Eq. 
hypergeomctric function ^: 
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= — Im 



The integrals in the last equation can be expressed in terms of the degenerate 
dte~ tp 
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where we used also the relationship 521 
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VP* ( 1, 



(B6) 



*(<z, b; z) = z 1 - b *(a - b + 1,2 - b; z). 

One can readily see that the integrals J% , Ji are monotonic functions of their arguments and have the following 
asymptotic expansions 
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Thus for regular part of the potential we obtain 
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while the oscillating part takes the form l|t).4|) written in Sec. IVII Eqs. I|B8|) and (|6.4|l were also derived in Ref. [24| 
using a different approach [53j. The sum in Eq. (|B8|I can be evaluated through the generalized zeta function in the 
following way 
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which finally results in Eq. I|tj.3f) . 

The oscillating part f2 osc (/x) of the thermodynamic potential can be also represented in terms of the generalized 
zeta function. For that one should use T(n + a) = J dss n+a ~ 1 e~ s and perform summation over n by means of the 
formula 
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Using now the formula 52] 



t u-x e -»x dx 

1 _ e -Px 



3" 



r(i/)c u, 



Re/i > 0, Ret/ > 0, 
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we find 
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This expression involves a dependence on the nonanalytical fractional part mod 
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In terms of the integer part 



2eB 



the total thermodynamic potential (a regular plus oscillating parts) can be written as |2£ 



~ , v NeB , 
SloiM) = - A M 
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(B14) 



The last two terms are in fact cancelled by the zero chemical potential (vacuum) part of the thermodynamic potential 

Eq. (£Q (for fj, > A). 

In fact, the expression Eq. (|B14|) can be obtained much easier if one uses Eq. I|4.2|l and writes Eq. I|6.4|l as 



NeB 
2tt 
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A + E - M « 
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Using now the formula for the generalized zeta function 



k-l 



C(z, v + k) = C(z, v) — J^(n + v) 
n=0 



(B16) 



one immediately arrives at Eq. I|B14(1 . However, the oscillating property of Ho(/i) is not so transparent in Eq. (|B14|) 
as it is in Eqs. (|(j.4fl and l|tj.6[l . On the other hand, the expression for f2 (A 4 ) through the generalized zeta-function is 
convenient for studying large field behavior. Indeed, for eB >> A 2 , /i 2 — A 2 we find from Ea. ljB14|l that in this limit 



Qo(m) 



NeB 
2tt 



%-A)(/x-A), 



(B17) 



which is precisely the contribution due to the lowest Landau level which is the only one to survive in the high field 
limit. Also, as is seen from Eq . (|B14jl . the oscillating term (first zeta-function in curve brackets) is necessary to take 
into account even in the high field regime in order to cancel (ei?) 3 / 2 term and obtain the correct linear term for Qq (/•*)• 
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APPENDIX C: CALCULATION OF n osc (» IN THE PRESENCE OF SCATTERING 

It turns out that due to the bad convergence of the integrals like (|5.7() and i|5.17[l with the Lorentzian distribution 
Pt(oj) this calculation can be done more easily if we directly substitute the DOS l|4.20|l in Eq. (|6.1|) and generalize 
the calculation made in Sec. I VII Choosing again for definiteness /i > 0, for the oscillating part of the thermodynamic 
potential we have 



a, 



*,r) = -^-A*-r*)i;-l / 
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exp 



(-2wky/j{x + u)) , 



where the variables w, u and 7 are defined in Eq. (|8.6|l . Using the representation 



(CI) 



2^ 7 t 3 / 2 
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we can perform the integration over x to get 
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dx sin(7rfcx) 
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(C3) 



where the variable v = fx 2 /eB is the same as in Eq. (16.51) . Hence the oscillating part of the thermodynamic potential 
can be written (omitting the monotonic term related to the first term in square brackets of the last expression) in the 
final form l|8.5|l using the functions of two arguments Ji(p,r) and ^{p, t) 



Vt(t 2 + 1) op 



dtVie-tP-T 



(C4) 



As is seen, the functions Ji(p, r), f) are monotonic functions of their variables. 

For large values of the parameter p we can obtain asymptotic expansion of the functions J\{p, r), J2{p, r) changing 
the variable t —* t/p and expanding then the denominator in powers of p, we get 



(-1)" f die'*'* 
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p2n+l/2 J f--2n+l/ 
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n+1/4 
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where we used the formula for the representation of the Macdonald function 
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and the relation K_ u (z) — K v (z), Since for half integer values of the index v 

(n + k)\ 



2: kl(n- k)\{2z) k 



(see Eq.( 



in |52j\ we obtain for the first few terms in the expansion 
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Similarly, for J 2 we get 

oo 
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FIG. 1: The magnetization M (in K/(Tesla cm 2 ) calculated for iV = 1) as a function of inverse field, B 1 , for three different 
values of chemical potential and A = T = T = 0. We use eB -» 200K 2 • B[Tesla]. 




FIG. 2: The magnetization M (in K/(Tesla cm 2 ) calculated for iV = 1) as a function of inverse field, B , for three different 
values of the gap A and for T = T = 0. We use eB — > 10 4 K 2 ■ _B[Tesla] and /i = 200K that are typical for graphite. 
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FIG. 3: The chemical potential (^/ef) 2 as a function of inverse field, B 1 , for three different values of the Fermi energy sf, 
F = 0.5K and T = A = 0. We use eB -» 10 4 K 2 ■ B[Tesla]. 
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FIG. 4: The magnetization M (in K/(Tesla cm 2 ) calculated for iV = 1) as a function of inverse field, B , for three different 
values of the chemical potential. We use eB -> 10 4 K 2 • S[Tesla], T = 0.5K, T = 5K and A = 0. 
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FIG. 5: The magnetization M (in K/(Tesla cm 2 ) calculated for N = 1) as a function of inverse field, B 1 , for three different 
values of the chemical potential. We use eB -* 10 4 K 2 • fl[Tesla], T = 0.5K, T = 15K and A = 0. 



